setwd("E:/2_5hmc_yjp_bam/ASM/varscan")
sel1=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=c("X2B_X1T","M8_M7","M6_M5")
sel1=c("M2_M1","M48_M47","M50_M49")
sel1=c("M18_M17","M20_M19","M22_M21","M40_M39")

library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)
library(progress)

sel1=paste(sel1,".snp",sep="")
snpdiff=c()
bayes_p=function(ref,var){
  df=data.frame(y=var,Ntrials=ref+var)
  null <- logit(median(df[,1]/df[,2]))
  formula = y ~ 1
  m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
  m=m1$marginals.fixed[[1]]
  lower_p1=inla.pmarginal(null, m)#lower_p
  upper_p1=1 - inla.pmarginal(null, m)#upper_p
  df=data.frame(y=ref,Ntrials=ref+var)
  formula = y ~ 1
  m2 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
  mb=m2$marginals.fixed[[1]]
  lower_p2=inla.pmarginal(null, mb)#lower_p
  upper_p2=1 - inla.pmarginal(null, mb)#upper_p
  
  bayes_pvalue=2 * (min(lower_p1, upper_p1,lower_p2, upper_p2))
  	###此处代码有点问题，因为我只想得到alt allele的beta0值，因此只取m1那部分就好了，修改后的代码同目录下名字为beta0_M1.R
  if(bayes_pvalue==2*(min(lower_p1, upper_p1))){
  coef <- m1$summary.fixed
  beta0=coef$mean
  lower_pvalue=coef$`0.025quant`
  upper_pvalue=coef$`0.975quant`
  }
  if(bayes_pvalue==2*(min(lower_p2, upper_p2))){
  coef <- m2$summary.fixed
  beta0=coef$mean
  lower_pvalue=coef$`0.025quant`
  upper_pvalue=coef$`0.975quant`
  }
  rt_num=paste(bayes_pvalue,beta0,lower_pvalue,upper_pvalue,sep="_")
  return(rt_num)
  
}

###原版本：
bayes_p=function(ref,var){
    df=data.frame(y=var,Ntrials=ref+var)
    null <- logit(median(df[,1]/df[,2]))
    formula = y ~ 1
    m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
    m=m1$marginals.fixed[[1]]
    lower_p1=inla.pmarginal(null, m)#lower_p
    upper_p1=1 - inla.pmarginal(null, m)#upper_p
    df=data.frame(y=ref,Ntrials=ref+var)
    formula = y ~ 1
    m2 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
    ms=m2$marginals.fixed[[1]]
    lower_p2=inla.pmarginal(null, ms)#lower_p
    upper_p2=1 - inla.pmarginal(null, ms)#upper_p
    return(2 * (min(lower_p1, upper_p1,lower_p2, upper_p2)))
}


for (i in 1:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
result=matrix(,nrow(file),ncol=2)
if(dim(file)[1]!=0){
fn="normal_tumor"
for(k in 1:2){
files=file[,c(paste(unlist(strsplit(fn,"_"))[k],"_reads1",sep=""),paste(unlist(strsplit(fn,"_"))[k],"_reads2",sep=""))]
for(j in 1:nrow(files)){
reads=as.numeric(as.matrix(files)[j,])
ref=reads[1]
var=reads[2]
	result[j,k]=bayes_p(ref,var)
	}
}

file$normal_rt=result[,1]
file$tumor_rt=result[,2]
file_name=paste0("./bayes_p_errobar/",gsub(".snp","",sel1[i]),"_bayes_errobar.txt")
write.table(file,file_name,quote=F,row.names=F,sep="\t")
}
}


##############################统计的基本信息用于填文章用的表
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p")
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel1=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=c("M18_M17","M20_M19","M22_M21","M40_M39")
sel1=paste0(sel1,".snp.bayes_p.txt")

col_names=c("FileName","Total_reads_up_10","Sig")
result=data.frame(matrix(NA,1,ncol = 3))
names(result)=col_names
result=result[-1,]
a3=c()
b3=c()
for(i in 1:length(sel1)){
  result_tmp=data.frame(matrix(NA,1,ncol = 3))
  names(result_tmp)=col_names
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="TRUE",]
file$unitID=paste(file[,1],":",file[,2],sep="")
file$total_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
file=file[file$total_reads>=10,]
a3=unique(c(a3,file$unitID))
result_tmp[,2]=dim(file)[1]
file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
b3=unique(c(b3,file$unitID))
result_tmp[,3]=dim(file)[1]
result_tmp[,1]=sel1[i]
result=rbind(result,result_tmp)
}


###以下为合并文件
rt=data.frame(as.matrix(snpdiff,1,length(snpdiff)))
names(rt)=c("unitID")
rt=tidyr::separate(rt,unitID,into=c("chr","pos","ref","var"),sep="_")
rt$unitID=paste(rt$chr,rt$pos,sep=":")
positions=as.character(unlist(rt$unitID))

i=1
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final=final[complete.cases(final$unitID),]
final$tumor_var_freq=as.numeric(gsub("%","",final$tumor_var_freq))/100
final$normal_var_freq=as.numeric(gsub("%","",final$normal_var_freq))/100
final$total_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)+as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
sels=which(final$total_reads>=10)
final=final[sels,]
final=final[,-8]

colnames(final)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")
}

for (i in 2:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final1=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final1=final1[complete.cases(final1$unitID),]
final1$tumor_var_freq=as.numeric(gsub("%","",final1$tumor_var_freq))/100
final1$normal_var_freq=as.numeric(gsub("%","",final1$normal_var_freq))/100
final1$total_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)+as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
sels=which(final1$total_reads>=10)
final1=final1[sels,]
final1=final1[,-8]
colnames(final1)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")

final=merge(final,final1,by="unitID",all=T)
}
}

###########读取文件并进行INLA分析

######ASE function要求结果file文件为count数，且顺序为
######twin1_control_ref,twin1_case_ref,twin2_control_ref,twin2_case_ref,
######twin1_control_var,twin1_case_var,twin2_control_var,twin2_case_var.

sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
#sel1=c("X2B_X1T","M8_M7","M50_M49","M52_M51");sel2=c("M28_M27")
#sel2=c("M30_M29","M26_M25","M35_M36")
#sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M42_M41","M44_M43","M48_M47","M50_M49","M52_M51","M12_M11");sel2=c("M28_M27","M30_M29","M26_M25","M36_M35")  ####包括抑郁症
sel3=c("M18_M17","M20_M19","M22_M21","M40_M39")

filekp=read.table("../bayes/AShM.txt",head=T,sep = "\t")
ref=paste(unlist(strsplit(sel1,"_")),"_reads1",sep="")
var=paste(unlist(strsplit(sel1,"_")),"_reads2",sep="")
file=filekp[,c(ref,var)]

ASE=function(x)
{
if (length(which(!is.na(x)))/4>=1){
sel=which(!is.na(x))
rct=x[sel]
ref=rct[1:(length(rct)/2)]
var=rct[(length(rct)/2+1):(length(rct))]
if (length(sel)/4>1){
df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4),x2=rep(c(0:(length(sel)/4-1)),each=2))
#header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1","x2")
#colnames(df)=header
formula = y ~ 1 + f(x2, model = "iid") + x1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
formula = y ~ 1 + f(x2, model = "iid")
m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
b10 <- exp(m1$mlik[2] - m0$mlik[2])
}else{
df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4))
#header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1")
#colnames(df)=header
formula = y ~ 1 + x1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
formula = y ~ 1
m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
b10 <- exp(m1$mlik[2] - m0$mlik[2])
}
return(b10)
}
}

Sys.time()
rt=apply(file,1,ASE)
Sys.time()
result=as.numeric(as.character(rt))

filekp$Bayes_in_DC=result
